Dynamical formation of domains in net-baryon density at the QCD phase transition 
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We consider the (3+1) dimensional expansion and cooling of the chirally-restored and deconfined 
matter at finite net-baryon densities as created in heavy-ion collisions. In our approach, chiral fields 
and the Polyakov loop propagate explicitly within a medium represented by a quark-antiquark fluid. 
The interaction between the fields and the fiuid leads to dissipation and noise, which in turn afi^ect 
the field fiuctuations. For the first time in a realistic dynamical study, we demonstrate the possibility 
of domain formation in a first-order phase transition of QCD matter. 
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The investigation of phase transitions in strongly- 
interacting matter during a heavy-ion collision is neces- 
sarily accompanied by nonequilibrium phenomena due to 
the fast collective expansion. At a first-order phase tran- 
sition, which is suggested by effective rnodels of QCD at 
large baryochemical potentials ^b [L, 2r], the rapid dy- 
namics should lead to strong supercooling. If nucleation 
times are large, a fast cooling will drive the system into 
the spinodal region, where the metastable state will de- 
cay via spinodal decomposition. This can lead to the 
formation of domains of the decaying phase embedded in 
the low-density background [J-Q. Since at small bary- 
ochemical potentials lattice QCD calculations predict a 
crossover transition [T], |8| , one may expect the first-order 
phase transition only at high values of fj,B starting from 
a critical end point (CEP). 

As the regime of high net-baryon densities will be ac- 
cessible in future experiments at FAIR [9] and NICA [lOJ , 
we will have the opportunity to directly investigate the 
first-order phase transition by looking for potential sig- 
nals stemming from the spinodal dynamics. 

In this letter we present a fully dynamical approach to 
the formation of domains at the first-order phase tran- 
sition during a heavy-ion collision. For this purpose 
we apply the recently developed model of nonequilib- 
riuETL (Polyakov-) chiral fiuid dynamics (N^FD, NP^FD) 
|llMl4l | , which extends earlier versions of chiral fiuid dy- 
namics [la ll6|. Here, the dynamics of the order pa- 
rameter fields and their respective fiuctuations are de- 
scribed by a stochastic Langevin equation coupled to 
the fluid dynamical expansion of the bulk matter. We 
adopt the Polyakov-Quark-Meson (PQM) model, which 
exhibits the expected thermodynamic phase structure of 
QCD at mean-field level [2] and beyond [17|. 

In the PQM model the mesonic fields a and 7? couple 
to the quarks and antiquarks with a coupling strength g. 



The temporal component of the color gauge field also cou- 
ples to the quarks, generating the Polyakov loop £ in the 
mean-field effective Potential K;ff which is represented as 



V,s = U{a) + U{1, T) + r!qq(a, ^, T, ^) 



(1) 



Here, U is the classical Mexican hat potential of the 
mesonic fields, W is a temperature-dependent effective 
potential of the Polyakov loop [2|, [l8|, |l9[ and the mean- 
field thermodynamic potential of the quarks reads 



n^^ = -2NfT , 



d^p 



(27r) 
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We used the abbreviating notation 






i(-) 



with the number of fiavors Nf — 2 and E^'^' — Eq — ^i, 
E^' — Eq + fjL, where /x = Ii^b/'^ is the quark chemi- 
cal potential . Throughout this work, we neglect pion 
fluctuations and concentrate on the sigma field and the 
Polyakov loop, which are the order parameters for the 
chiral and deconfinement transition. The energy of the 
quarks Eq — Jp^ -\- mi depends on the sigma field via 

the dynamically generated effective mass ruq = g\a\. In 
the mean-field approximation, the chiral and deconfine- 
ment transition coincide for all values of fx. We also find 
a common CEP at (Tcp, ^cp) = (152, 160) MeV. 

For effective models of the present type it was found in 
20l I2I I that nucleation times are much longer than typ- 
ical expansion times of the fireball, which implies that 
the system will strongly supercool. This conclusion was 



confirmed within the present model in 1^ |lj| , thus giv- 
ing us the reason to expect observable signals stemming 
from spinodal decomposition. 

We now briefly describe the nonequilibrium Polyakov- 
chiral fluid dynamics model (NP^FD). In this approach 
the sigma field is propagated according to a Langevin 
equation 



d^d^'a 



VcrdtO- 



Sa 



(5) 



which besides the standard mean-field contribution con- 
tains the damping term with friction coefficient r]^ and 
the stochastic noise field ^a- For m„{T) > 2mq{T) = 
2gacq{T) the damping coefficient due to the decay of 
the zero mode sigma field into a quark-antiquark pair 
(ct — > q + g) is given by [11| 
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(6) 



Since we are mainly interested in the long-range fluctu- 
ations of the sigma field, we use the damping coefficient 
([6]) for all modes in the Langevin equation. At lower 
temperatures it vanishes because this decay is kinemat- 
ically forbidden as ma{T) < 2mq{T) — 2gacq{T). How- 
ever, the sigma field is further d amp ed due to interaction 
with the hard chiral modes [221 123| , which can be repre- 
sented by the reaction a ■<-> 27r. To account for them we 
use the value 77 = 2.2/fm [24] for the kinematic range 
2mq > ma(T) > 2-171^^. 

The stochastic field in the Langevin equation ([5]) is 
given by a Gaussian distribution with zero mean 



(C^(i))c = 0. 



(7) 



The amplitude of the noise is determined by the 
dissipation-fluctuation theorem 



2t) 



'^' (8) 

The dynamics of the Polyakov loop t is described by 
the relaxation equation 






ii 



(9) 



with a constant damping coefficient rji — 5/fm. This 
value can only be estimated very roughly, but the results 
are largely independent of its specific choice [1J|. The 
stochastic noise ^i is assumed to be Gaussian, too, and 
we impose the dissipation- fiuctuation relation 14 1 



iUt, x)Ut', S'))^T^ - ^^t - t')5{^ - x')2T^tT . (10) 

In our dynamical model, the expansion of the fireball 
is described according to energy-momentum and baryon 
number conservation equations 






(11) 
(12) 
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FIG. 1. Trajectories of the hydrodynamic simulation in the 
phase diagram for different initial conditions. Nonequilibrium 
effects lead to an overshooting of the first-order phase transi- 
tion. The dashed gray lines denote the spinodal region. 



where the source term S'^ accounts for the energy- 
momentum exchange between the fiuid and the fields. 
It is dominated by the energy dissipation from the sigma 
field and Polyakov loop to the fiuid, controlled by the 
damping terms in Eqs. ([5]) and ^. This energy trans- 
fer is largest during the relaxation from the supercooled 
state. 

Since the fields a and t are not necessarily at their 
respective equilibrium values, the equation of state will 
depend explicitly on their local values cr(x), t{x). This 
is different to standard fluid dynamical simulations using 
equilibrium equations of state. The corresponding equa- 
tion of state of the fluid is obtained from the generalized 
thermodynamic relations 



n{a,e,T,fi) = — . 



(13) 
(14) 

(15) 



where a and £ are functions of time and coordinates due 
to their own dynamics. 

In the numerical implementation we solve the coupled 
equations ([5]), © and the (3+1) dimensional fluid dy- 
namic equations ([TTj). (IT^ . as it is described in detail 



13j, |lj]. As initial states we consider spherical flre- 



balls with various initial values of T and /x to probe the 
crossover, CEP and first-order transitions. The order pa- 
rameter fields are initialized in local thermal equilibrium 
for the initial T-/L(-profiles. 

The noise fields are sampled for each new time step 
according to the Gaussian distribution, given by Eqs. (|8]) 
and (fTU]) . For each grid cell (i G 1 . . . N) we average the 
initially uncorrelated noise variables over a surrounding 
volume V^°" = 1/to^ and V^°" = l/ra\ and attribute 
this average value to the noise in the respective grid cell 
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FIG. 2. (Color online) Relative net-baryon number density at z = for a first-order transition at t = 6 f m (a) and t = f 2 fm 
(b) and for a transition through the CEP for i = 6 fm (c) and i = 12 fm (d) One can see that droplets of high density are 



formed at the first-order phase transition, while the density fluctuations remain small during an evolution through the CEP. 



i. Thus, the noise values in neighboring cells i and i + 1, 
which share large parts of the averaging volume, are cor- 
related. This averaging reduces the variance of the noise 
and we need to recover the original local variance by 
rescaling. By this procedure we avoid an artificial corre- 
lation length equal to the cell spacing. In 1J| we have 



found that during the simulation, the correlation length 
reaches values 1.5 — 2.0 fm around the CEP and remains 
small, at about 0.2 — 0.4 fm, in the case of a first-order 
transition. 

The trajectories of the dynamic simulations in the 
T-^-plane are presented in Fig. [T] together with the cor- 
responding phase boundary and the spinodal lines. The 
averages (T) and (/z) are calculated over a central box of 
1 fm . The bending of the curves towards larger values of 
/i is due to the rapid increase in the quark mass at the chi- 
ral transition. Similar behavior is already inherent in the 



equilibrium isentro pes , which follow the phase boundary 
towards larger n [ij, |25| . In our dynamical simulations we 
see that at larger chemical potentials the bending only 
occurs after the trajectories have crossed the lower spin- 
odal line. This behavior has a simple explanation: After 
crossing the first-order phase transition line, the system 
is still not able to overcome the potential barrier but gets 
trapped in a metastable chirally restored and deconfined 
state. After that barrier has vanished at the lower spin- 
odal line, the system can roll down to the global min- 
imum where the constituent quarks acquire their mass. 
As seen in Fig. [1] at highest fi, one can observe a re- 
heating of the quark fluid, an effect that we have alr eady 
found in calculations at zero chemical potential 13|, [ij] . 
In summary we see that due to nonequilibrium effects 
the system supercools and spends an extended amount 
of time in the spinodal region. This may facilitate the 
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FIG. 3. (Color online) Azimu thal distribution of the baryon 
number density after t — Q fm|(a)|and f = 12 fm (b) for sev- 



eral transition scenarios. We see that strong inhoniogeneities 
develop at the first-order phase transition. 



process of spinodal decomposition and therefore the de- 
velopment of interesting signals of the first-order phase 
transition. 

In Fig. [21 we show two-dimensional maps (through the 
center of the fireball) of the relative net-baryon num- 
ber density n/{n), for evolutions through both the first- 
order phase transition and the CEP. They correspond to 
the most right solid line and the dashed line in Fig. [1] 
Here, (n) is the volume average over all cells with n > 0. 
We see that during the expansion the fluid evolves inho- 
mogeneously in the case of a first-order transition, cre- 
ating droplets of high baryon density at the periphery 
of the system. This effect is not observed at the CEP 
or crossover transition, where the spherical symmetry is 
preserved and no strong fluctuations occur. 

It is clear that such inhomogeneities should lead to 
large fluctuations in the angular distribution of the net- 
baryon number dN/dcfi, where (j) is the azimuthal angle 
around the z-axis. Such distributions are shown in Fig. 
[3] for i = 6 fm and i = 12 fm. Here one can clearly see 
strong fluctuations of the net-baryon number as a func- 
tion of 4> at the flrst-order transition in contrast to the 
CEP and crossover. One can notice the correlation of 



bumps and deeps between the plots in Fig. [3] that in- 
dicates that the domains preserve their identity during 
this time interval. From Fig. [2] one can see that individ- 
ual domains move radially outwards while their relative 
strength is amplified, since the system spends a long time 
in the spinodal region, see the corresponding trajectory 
in Fig. [T] From these figures we expect that after freeze- 
out the particles stemming from the domains will show a 
characteristic azimuthal spectrum leading to an enhance- 
ment of the higher fiow harmonics. 

In conclusion, we have demonstrated for the first time 
in a realistic dynamical study the formation of domains 
in net-baryon density at the first-order phase transition. 
These domains do not result from initial inhomogeneities 
but are due to the non-deterministic evolution of the sys- 
tem during spinodal decomposition. Within the setup of 
NPxFD we have shown that a clear signal for a first- 
order phase transition can develop and survive during a 
realistic simulation of the expansion dynamics which is 
imprinted in net-baryon azimuthal fiuctuations. These 
non-statistical fluctuations should be observed in single 
events. For a set of many events the prominent signa- 
ture would be an enhancement of higher flow harmonics. 
By varying the beam energy one can get closer to the 
CEP where the signals of the first-order phase transition 
become weaker, until they finally cease. Being able to 
perform an energy scan, which covers small enough re- 
gions of iiB, would thus offer the possibility to discover 
and locate the CEP by looking where the signals of the 
first-order phase transition vanish. 
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